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Abstract. Vortices of several condensed matter systems are predicted to have zero- 
energy core excitations which are Majorana fermions. These exotic quasi-particles are 
neutral, massless, and expected to have non-Abelian statistics. Furthermore, they 
make the ground state of the system highly degenerate. For a large density of vortices, 
an Abrikosov lattice is formed, and tunneling of Majorana fermions between vortices 
removes the energy degeneracy. In particular the spectrum of Majorana fermions in 
a triangular lattice is gapped, and the Hamiltonian which describes such a system is 
antisymmetric under time-reversal. We consider Majorana fermions on a disordered 
triangular lattice. We find that even for very weak disorder in the location of the 
vortices localized sub-gap modes appear. As the disorder becomes strong, a percolation 
phase transition takes place, and the gap is fully closed by extended states. The 
mechanism that underlies these phenomena is domain walls between two time-reversed 
phases, which are created by flipping the sign of the tunneling matrix elements. The 
density of states in the disordered lattice seems to diverge at zero energy. 
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1. Introduction 

Exotic states of matter are among the most intriguing topics in the field of condensed 
matter physics. One class of these exotic states are the two-dimensional (2D) systems 
in which quasi-particles follow non-Abelian quantum statistics [T] . The search for such 
systems is driven both by their unique properties and by their potential application to 
topological quantum computing [2]. 

In a non-Abelian state the ground state is degenerate when quasi-particles are present, 
and the degeneracy increases exponentially with the number of quasi-particles. Perhaps 
the simplest way of obtaining such a degeneracy is by having quasi-particles that carry 
Majorana fermionic excitations [3]. Majorana fermions (MFs) are expected to appear 
as zero-energy excitations in the cores of vortices in a layered p x + ip y superconductor 
- such as proposed for Sr 2 Ru04 pi], in 2D systems that can be mapped onto such a 
superconductor - the v — 5/2 fractional quantum Hall state [5], on the surface of a 
topological insulator that is in proximity to an s-wave superconductor and an insulating 
ferromagnet [6], and at hybrid structures of semiconductors and superconductors 
[7]. Furthermore, they are expected to form in ID systems in proximity to s-wave 
superconductors [8] , and in several other systems (9j [10] ■ 

A variety of experiments have been proposed in order to probe the predicted MFs, 
based on the their unique properties. To mention a few examples: zero-energy 
excitations can be observed by performing STM measurements at the vortex core 
[TT] IT2] ; the degeneracy of the ground state affects the thermodynamical properties 
[13] HI] ; non-locality of an electron in a Majorana state has a signature in tunneling 
between two vortices [TS1 HE]; and interferometric experiments are able to probe the 
non-Abelian statistics [2] HT] HE] CEH] • 

The suggested thermodynamics and interferometry measurements are based on the 
unique many-body properties of the MFs, but assume that the MFs are localized at 
theirs positions. The MFs appear at vortex cores, which are expected to rearrange as 
an Abrikosov lattice at high enough density. The lattice order and the small tunneling 
amplitude between neighboring vortices remove the ground state degeneracy and form 
a band of low-energy excitations. On one hand, this band may serve to probe the 
existence of the MFs. On the other hand, it may conceal signals of suggested many- 
body measurements, especially controlled adiabatic processes, such as interferometry. 

Some previous works have analyzed clean periodic square, triangular [20J and 
honeycomb [21J lattices of MFs, and found the electronic conductivity associated with 
tunneling between vortex cores. Other works have considered some aspects of random 
square [10] and honeycomb (22] lattices. In this paper, we consider disordered triangular 
lattices. This lattice breaks time-reversal symmetry: it may be mapped onto a tight- 
binding model of electrons on the same lattice with each plaquette being pierced by 
±1/4 magnetic flux quantum. The sign of the flux is determined by the sign of the 
tight-binding coupling term. The spectrum of the perfect lattice is gapped. 

Our main finding is that disorder in the sign of the tunneling amplitude between 
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neighboring sites creates a peak in the density of states (DOS) at zero energy. We show 
that at weak disorder zero-energy localized sub-gap states are created, and coupling 
between these states creates a DOS close to zero energy. At strong disorder the sign 
flipping creates domain walls between regions of the two time-reversed phases, and chiral 
modes appear along these walls. When the signs are random, we find the domains to 
be narrow and the spectrum of excitations to be characterized by strong dependence on 
the geometry of the walls. Plausibly, that is the source of the zero-energy DOS peak in 
the highly disordered system. The transition between these two limits is a percolation 
phase transition, taking place when the probability of flipping a sign is around 0.15. We 
find the DOS to diverge at zero energy at intermediate and strong disorder. 

The paper is organized as follows: In section [2j we define the MFs lattice model, 
and show the numerical DOS for a disordered lattice. In section [3j we examine how 
low-energy excitations may emerge from localized modes in the limit of weak disorder 
and from extended interfaces in the limit of large p. Section [4] discusses the dependence 
of the excitation spectrum on the strength of the disorder. Section [5] summarizes the 
results and compares them to previous works. 

2. A disordered triangular lattice 

Our interest here is in MFs that form a triangular lattice. An MF is defined as a self- 
adjoint operator 7 = 7' that satisfies fermionic anti-commutation relation {7^ 7.,-} = 5y. 
A standard complex fermionic operator can be constructed from MFs by superposing 
an even number of them. For example ip n = (7, + i , jj)/\/2 for two MFs satisfies the 
standard anti-commutation relations {VvV'm} = $nm and {ipn,ipm} = 0. 

We consider MFs that are solutions of the Bogoliubov de-Gennes (BdG) equation 
in the presence of vortices in the superconductor. Each vortex in the superconductor 
carries a signle MF, whose wavefunction is exponentially localized around the vortex's 
core. Overlaps of the MF wavefunctions of neighboring vortices result in tunneling 
matrix elements. Vortices in superconductors and quasi-particles in clean quantum Hall 
systems are arranged on a lattice; thus the MFs are also arranged as a lattice, and their 
dynamics may naturally be described by the tight-binding model. 

The particle-hole symmetry of the BdG equation implies that a single MF is a zero- 
energy solution. Hence, in the tight-binding Hamiltonian the on-site energy of the MFs 
is zero, and there are only hopping terms tij. Moreover, for the Hamiltonian to be 
Hermitian, Ujji'jj = (tijjijj)^ = —t*jlijj, implying a purely imaginary hopping term. 
By assuming discrete symmetry of the lattice to translations, we can write the simple 
Hamiltonian 



where (ij) are nearest neighbors and s^- = — Sji = ±1. 

Any element Sij is gauge dependent, because the 7^ operators are defined up to an 
overall sign. However, the product of s^-'s along a path creating a closed loop is gauge 
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independent. It has been shown in [20] that for any lattice whose plaquette is a polygon 
of n vertices, the product of Sy around each plaquette corresponds to the plaquette 
enclosing n/4 — 1/2 flux quanta. Therefore the product of the hopping terms along the 
bonds that create the plaquette is —i n t n , and the product of the s^-'s along this path 
is —1 (the direction of the path is chosen to be aligned with the chirality of the order 
parameter) . 
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Figure 1. (a) One possible gauge of the tight-binding model of MFs on a triangular 
lattice, which is expressed in equation pj). An arrow from site i to site j means = 1. 
Note that the lattice is split into two sublattices A and B, and that the counterclockwise 
product of the Sj,-'s around each triangle equals to —1. (b) The corresponding gapped 
spectrum, with E„ ap = 1.73. 



In particular, in the triangular lattice each plaquette encloses a quarter of flux 
quantum, and the product of the hopping terms equals it 3 , revealing a time- reversal 
anti-symmetry. Moreover, the unit cell which encloses a flux quantum is a parallelogram 
of four neighboring triangles, leading to the lattice being composed of two sublattices. 
One possible gauge is illustrated in figure [l|a), which also depicts the two sublattices A 
and B. In this gauge 

H = it)] (-7A,i7Ai+77i + 7B,i7B,i+T7i + lB,ilA,i- m 
i 

+lA,ClB,i- m ~ "fB,aA,i+ Vl -r, 2 + nfA,ijB,i+ Vl -r, 2 ) ■ (2) 

Let us assume a clean periodic lattice of L\ x L 2 sites, with L 2 even. We can define the 
Fourier transform T a k = Yli e ~ lkXi 1a,i> where a = A,B and k = 27c(mi/ L\, m2/2L2) for 
rrii = 0, ...,Lj — 1 (i — 1,2) [20]. These transformed operators are complex fermions, 
which satisfy {r a fc , fc ,} = <5 fcfc /. Note, however, that for k = (0,0), (tc,0), (0,7r/2) and 
(7r,7r/2) they are MFs. By denoting spinor = (f^ k ,^ B k ) the Hamiltonian can be 
expressed as 

H = 2t f k (sin k 2 a x - cos(A;i - k 2 )o- y - sin ki<j z ) f fc , (3) 

k 

where the Pauli matrices act on the sublattice space. The resulting spectrum has an 
energy gap of 2tE gap , where E gaiP = 1.73, as depicted in figure [T](b). 
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The hopping elements between neighboring MFs are very sensitive to the inter-vortex 
separation. For example it was shown that for a p x + ip y superconductor |23j : 



t ~ cos 



(k F r + 



-r/t 



(4) 



where kp is the Fermi momentum, r is the inter-vortex distance, and £ is the coherence 
length, which is usually larger than kp. The exponential decay and the oscillations 
were found to occur also for the v = 5/2 case [21], and are likely to appear in all 
realizations. Therefore small deformations of the Abrikosov lattice of order k-^T 1 will 
produce fluctuations in both the amplitude and sign of the hopping terms ty. 
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Figure 2. The disorder-averaged DOS (D) as a function of the energy E of a periodic 
30 x 30 lattice with (a) uniform and (b) random hopping terms. In systems with 
hopping terms that are random both in sign and in magnitude (blue) and in systems 
where only the signs are random (green) there is a sharp peak at zero energy, while 
random amplitudes (red) only smear the spectrum of the clean system, but do not 
close the gap. The DOS is averaged over 10,000 realizations. 



We find numerically that random hopping terms ty = tsy, where s^- is uniformly 
distributed in the interval [—1, 1], close the energy gap of the uniform system. The DOS, 
which appears in figure |2j shows that not only the gap is closed, but a peak emerges at 
zero energy. We can distinguish between random amplitudes and random signs of the 
hopping terms, i.e. |sy| ~ U[0, 1] or s^- = ±1 in equal probability, respectively. The 
DOS of these two cases, which are also depicted in figure [2j clearly shows that the zero- 
energy peak in the DOS appears only due to random signs, while random amplitudes 
merely smear the spectrum without closing the energy gap. We are mostly interested 
in the zero-energy peak of the DOS, and therefore in the following we will focus on the 
case where disorder appears in the sign of the hopping terms. 

The sharp peak that we found numerically for the density of states of finite systems 
close to zero-energy naturally raises the question of the way the DOS scales with the 
size of the system in the thermodynamic limit. Increasing the system size L obviously 
increases the DOS. If the lowest energies decay faster than 1/L 2 , then the zero-energy 
DOS will increase faster than L 2 , and the zero-energy DOS per unit area would diverge 
in the thermodynamic limit. Figure [3|a) exhibits the energy of the three lowest states 
Ei,E 2 and E 3 , averaged over disorder, as a function of the system size. In the given 
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Figure 3. (a) The dependence of the disorder-averaged energies of the three lowest 
states E\,E2 and E3, as a function of the system size L 2 , on a log- log scale. The 
energies (E n ) were averaged over 1000 realizations. Linear fitting gives (E n ) oc L~ 21S 
for all three energies, (b) The disorder-averaged DOS (D(E)) with first (t) and second 
(i 2 ) nearest neighbors hopping, both with random signs, for t 2 /t = 0,0.5,1. As £2 
approaches t the DOS approaches a semicircle. The DOS is of a 30 x 30 lattice, and 
was averaged over 1000 realizations. 



It is instructive to compare the DOS we find for the case of random nearest-neighbor 
hopping to that we find for the case of a random Hermitian matrix of imaginary terms. 
The DOS of the latter is a version of a semicircle distribution, with a Delta peak at 
zero energy, due to the symmetry of the spectrum with respect to reflection about 
zero energy [25]. When we add second nearest-neighbor hopping with random signs to 
the Hamiltonian, the DOS indeed gets closer to the semicircle. Figure |3](b) shows the 
spectrum for several ratios of t 2 , the amplitude of the next-nearest-neighbor terms, to t. 
We note, however, that we numerically find the gap to close even for the clean system 
when t 2 /t « 0.58. 

Having established that the DOS of the clean lattice is characterized by an energy gap 
limited by two singular peaks, and the lattice where the sign of the hopping terms is 
random is characterized by a zero energy peak above a gapless background, we examine 
the evolution of the DOS with p, the probability of flipping the sign of a hopping 
term, which increases from (the clean lattice) to 0.5 (the random sign limit). Figure 
|4]^a) shows how increasing p makes the singular peaks of the clean system spread, and 
gradually creates sub-gap states. 

In the next section we study the way the gap is closed, and especially the appearance 
of the zero-energy peak. We address these questions by an examination of two limits. 
In the weak disorder limit, we present the minimal disorder configuration that results 
in a zero-energy state. In the strong disorder limit, we examine the low-energy states 
that are formed along domain walls between regions in which the signs of flux piercing 
the plaquettes are opposite. 
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Figure 4. The disorder-averaged DOS (D(E)) of a periodic 70 x 70 lattice as a 
function of the probability p for flipping the sign of a hopping term. For p = the 
DOS is that of figure [2]ja). Two small peaks at E « ±1, that grow linearly with p, 
belong to localized states which are created around a single isolated flipped hopping 
term. The zero-energy peak grows slowly, but survive at the disordered system. The 
DOS is averaged over 100 realizations. 



3. The origin of the zero-energy peak in the density of states 

We have seen that the clean triangular lattice of MFs is anti-symmetric under time 
reversal, and has an energy gap. Depending on the sign of the hopping tunneling 
elements, the MFs may be mapped onto electrons in a triangular lattice pierced by 1/4 
or —1/4 flux quanta per plaquette. These two cases are characterized by two opposite 
Chern numbers: by writing the k th component of equation ([3j as = 2th & ■ cr, it is 
easy to see that the Majorana band is characterized by a non-trivial Chern number 

• M f d2k 1 (u dh >° dh * 
v = srgn(t) / — — - h k ■ x 



47r \hh\ 3 \ dki dk 2 

2 sin 2 (A;i — k 2 ) + cos k\ cos k 2 cos(/ci — k 2 ) 



dfci / dk 2 - 

■7T J-7r/2 



-tt/2 Att [sin 2 ki + sin 2 k 2 + cos 2 (k 1 - k 2 )] 3/2 
= sign(t). (5) 

The disorder-induced reversal of signs of hopping terms may revert the sign of the flux 
in some of the plaquettes in the lattice. The systems would then have islands of one 
phase separated by lines of interface from the bulk of the other phase. An infinite 
interface between two phases of different Chern numbers is accompanied by gapless 
modes j26j ETJ EE] . For finite islands, one may naively expect finite-size quantization to 
induce a gap in the energy spectrum, and thus low-energy excitations to require large 
islands. This expectation is only partially valid. As we explain below, for the Majorana 
lattice, large islands are associated with low-energy excitations at their edges, but small 
islands may carry low-energy and zero-energy modes as well. 

When the probability p of flipping the sign of a hopping term is very small the flipped 
terms are dilute, and most of them are isolated. A two-triangles island, which is created 
by flipping the sign of a single hopping term, results in two localized sub-gap states with 
energies of approximately ±1.1. Theses states are the source of the peaks in the DOS 
at ±E rj ±1.1, and indeed in the limit of small probability the DOS associated with 
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these peaks {D(E ~ ±1.1)} oc p. Such islands do not, however, contribute to the DOS 
close to zero energy. 

As we show in Appendix A , flipping three hopping terms with a common vertex 
creates two localized states with energy that is either zero or exponentially small with 
the system size. We also show that this is the minimal way of creating zero modes. 
Thus, for small p the zero-energy DOS is dominated by the probability of creating such 
configurations, and should scale like p 3 . The two modes are, however, split when there 
is a single flipped hopping term at a distance r from that vertex. The exact splitting 
depends on the orientations of the bonds with flipped signs, but we found numerically 
that it can be well approximated by ±E t (r) m ±6e _2r . For E t to be much smaller 
than the typical finite-size energy splitting, that scales as L~ 2 , we need r ^> logL. A 
given radius r encloses approximately 3irr 2 hopping terms around the vertex, and the 
probability that all these terms are unflipped is (1 — p) 3nr2 . Thus, the p 3 dependence 
of the zero-energy DOS is limited to small values of p <C log -2 L. For example, for a 
70 x 70 lattice we get p < 0.01. 

As p gets large, the plaquettes in which the flux is reversed connect to one another, 
and the system is split into domains with opposite fluxes. To understand the excitations 
spectrum in this limit, we first examine the spectrum of excitations associated with the 
various possible interfaces of two large domains of opposite fluxes. Then, we examine 
the statistical distribution of such interfaces as a function of the flipping probability p. 

At the interface between a macroscopic region of Chern number v — ±1 and the 
vacuum, a chiral gapless mode must appear [26j E3 EE]. In order to explicitly find this 
mode, we first replace r]\ and 772 in equation ^ by x and y for simplicity, and denote 
the lattice sites by x. We assume rotational symmetry along y, and define the operators 
r a ,fcr = Yl y e ~ lky Ja,x with a = A, B and k = nm/L y , where m = 0,...,L y — 1. Each 
r a ,fcr represents a wavefunction that is localized at x but extended at y. With these 
operators the Hamiltonian becomes 

H = ~ l E [ r U e ~% + ^) r M+i + rUe% - u7,)r M _! - r^(2sin^.)r M ] ,( 6 ) 

Given a right edge at x = 0, the edge eigenmode of (|6p for x < is 
E^cA'^M)^,^, where b k = ie" lfc tan(fc/2) w i\k < 1 and (a, b)T k>x = 
o^A,k,x + bT B)kx . For k — it reduces to (1, l)r . The dispersion of the edge 
mode is E k = 2sinfc, with E k w 2k — 2irm/L y at low energy. Note that the 
wavefunction is localized in the x direction with a localization length of the order of 
l/log(L y /7rm) ~ l/log-E, which depends weakly on the energy and the system size. 

A triangular lattice with one periodic coordinate can have either flat or zigzag edges, 
depending on the way the periodic boundary condition is taken, as illustrated in figure 
[5} The Hamiltonian ^ describes a flat edge. In a zigzag edge, we also get a chiral 
edge mode, with the dispersion E k ~ y/2k and b k ~ (3 — 2\/2) — (8 — ll/\^2)k 2 , which 
gives a localization length that is approximately independent of energy, at low energies. 
We can see that in spite of the differences, both edges show a similar behavior of linear 
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dispersion and exponentially localized wavefunction. Henceforth we will assume the 
periodicity of equation 




Figure 5. Two types of edges in the triangular lattice with one periodic coordinate: 
(a) flat and (b) zigzag. With the periodic coordinate of (a) there are two types 
of domain walls between P-phase (azure) and N-phase (pink) that preserve the 
translational symmetry: (c) flat and (d) sawtooth. 

An interface, or domain wall, between the two phases of opposite Chern numbers is 
expected to hold two chiral edge modes, similar to the single mode that exists along 
the interface between any of the two phases and the vacuum. We denote the unflipped 
phase by P and the flipped phase by N, representing their positive and negative fluxes 
respectively. A periodic lattice with an interface between a P-phase and an N-phase 
may have two kinds of domain walls that preserve the translational symmetry: flat and 
sawtooth, as depicted in figure |5j 

A flat domain wall is created when o~ y — > —o y in the hopping terms for x > 0. The 
two chiral modes of such a domain wall have the same dispersion and amplitudes bk as 
the flat edge mode, where in the N-phase the amplitudes are b k *, as expected from the 
time reversal. Only now one mode has non-vanishing amplitudes at x — 2n, while the 
other mode has non- vanishing amplitudes at x — 2n + 1. A sawtooth domain wall is 
a result of flipping the cr^'s. It shows essentially a similar behavior to the flat domain 
wall, but with Ek ~ 1.95(A; ± k z ), where k z = arccos (VE/2 - 1/2). A gain, we see that 
the geometry of the domain walls has only a small effect on the low-energy physics. 

If an 'island' of a roughly circular shape of the N-phase is surrounded by the P-phase 
lattice, then the two chiral modes will be localized along the perimeter of the island, 
and their dispersion will be E m ~ m/M, where M is the perimeter of the island. Such 
an island will close the energy gap only if it becomes macroscopic. 

In a disordered system where the sign of the hopping term Sy is random, islands of 
the N-phase will be created in various shapes and sizes. A way to parameterize the 
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morphology of an island is by the perimeter-area relation. In approximately rounded 
shapes M oc A 1 / 2 , where A is the area of the island, while in a narrow strip M oc A. 
Figure [6|a) depicts the distribution of the perimeters as a function of the areas of 
random islands in a 30 x 30 lattice, where the probability of flipping a sign is 1/2. 
The area of an island is the number of N-triangles that compose the island, where two 
N-triangles are defined to belong to the same island if they share at least a vertex. 
The perimeter is composed of the P-triangles which have a common vertex with the 
island. For islands with A > 10 the mean perimeter can be excellently approximated 
by (M(A)) rj 2.3AA + 15.6, which means that the islands look like snakes or gossamer. 




20 40 60 80 100 
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Figure 6. (a) The distribution of the perimeters M of the N-phase islands with a 
given area A, where the islands were created by hopping terms with random signs. The 
error bars are located at the mean values (M), with the width of a standard deviation. 
The mean value is excellently fitted to a linear curve (red line). The distribution 
was taken from 10,000 random realizations of a 30 x 30 lattice, (b) A typical islands 
configuration of N-triangles (pink) that are created by random flipping of hopping 
terms (red), with a complicated geometry of rhombi and triangles. In this realization 
0.14 of the hopping terms are flipped. 



The smallest island that may be created by flipping the sign of a single hopping term 
is composed of two triangles, and has a rhombus shape. The simplicity of the shapes 
is not preserved for larger islands. Flipping two neighboring hopping terms results in 
two flipped N-triangles, but those share only a vertex. Larger islands are typically of 
complicated shapes, composed of rhombi and triangles, as illustrated in figure [6^b). 
Complicated islands provide complicated dispersion and DOS, as we will demonstrate. 
In order to get a feeling to the resulting dispersions, we now focus on three types of 
narrow islands that preserve the translational symmetry, i.e. strips of N-triangles. 

The method we use for approximating the low-energy dispersion of the strips is first 
to find zero-energy modes with k = 0, and then extract the dominant /c-dependence 
employing first order perturbation theory. For the clean system the Hamiltonian ^ 
may be recast as 

H =2tJ2 T aA- a y + 1(7 z) T o,x+i, ( 7b ) 
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5H k = 2t sin \ J2 (e-«/ 2 rl x i* y r k>x+1 - e^ 2 T{ x ia y T k>x ^ + 2 cos ^T{ x a x T k/x )j . (7c) 

For any strip s, this Hamiltonian has to be modified to the Hamiltonian H s , which may 
be decomposed similarly to Hq and SH k . The k = zero modes of H s will be denoted 
by q = rf. They satisfy [H°,Tj} = 0. In the cases we will explicitly consider here, 
there are two zero modes per strip. The matrix elements of SH k between these zero 
modes give the leading order of the /c-dependence. Therefore after finding Tf and Tjj, 
we diagonalize the matrix [5Hk], where [8Hk]ij = [5Hk,T*]} 
The first and simplest strip is the flat strip, which is depicted in figure [7]^a), and is 



i^Q at has two zero modes: 



rr = _(i,i)r 0i0 , 



created by flipping the sign of the cr z 's in equation <\7b\j between x = and x — 1. Now 

ie 
1 

7! 

rjS rt = ^(i ) -i)r ,i. (8) 

The matrix elements of 5Hk in the zero modes space are 

[5H k ] aat = -t sin k (l, - tan ^, 2^ • r, (9) 

with r being the vector of Pauli matrices. The resulting corrections to the zero energies 
are 5E aa,t ±y5tk. These modes are essentially the domain wall modes, which are 
modified due to their hybridization, and similarly their DOS is a constant. 




Figure 7. The three examples of narrow islands of N-phase, whose chiral sub-gap 
modes show three kinds of dispersion: (a) a flat strip with oc k, (b) a 'zipper' strip 
with Ef. oc k 2 and (c) a sawtooth strip with E^ oc k . 



The second example is a 'zipper' strip, shown in figure [Tj^b) , which is a result of 
replacing o~ y of Equations (7b) and (7c) by io~ x for x — — 1 and by — \o x for x — 0. Two 
zero modes again appear 

rf p = -L [(i, i)r 0> - 2 + (o, 2)r , - (1, -i)r„, 2 ] , 



-\zip 



2^2 
1 



[(i,i)r ,-i + (i,-i)r ,i] 



(10) 



which yields [5H k } zip = -y/2t sin 2 |cr r The corrections to the energies are now quadratic 



5El ip ps ±^tfc 2 . Note that these modes are wider, and that their DOS diverge as E l / 2 . 
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The third example is the sawtooth strip, which is created by multiplying the hopping 
terms between x — — 1 and x = by o~ z , and is shown in figure [7](c). The zero modes 
are now 

rr = ^(i,-i)r ,o, 

rr = ^[(i,i)ro,-i + (i,-i)r 0l i]. (11) 

Unfortunately, [5Hk] saw = —t(l + r z ) sin k, which means that the first order perturbation 
theory gives the correction only to the first mode 5El aw w —2tk, while the second mode 
requires higher orders. An explicit solution of this mode gives 5£| aw w fc 3 ; thus this soft 
mode makes the DOS diverge as E~ 2 ^ 3 . 

We can conclude from this section that along edges and domain walls chiral modes 
appear with linear low-energy dispersion. These modes close the energy gap only 
for macroscopic domain walls, which will give a constant DOS. However, when the 
domains become narrow strips, the dispersions of the chiral modes strongly depend on 
the geometry of the strip, and their DOS tends to diverge at zero energy. And since 
most domains which are created by random hopping terms are composed of narrow 
strips, this may explain the zero-energy peak raising above the uniform background in 
the DOS. 



4. Percolation phase transition 

In the previous section we have seen that chiral modes which are localized along stripe- 
shaped domains may lead to a zero-energy peak in the DOS when the signs of the 
hopping terms are random. Such modes provide low-energy states (states whose energy 
scales inversely with the size of the system) only if the length of the domain is of the 
order of the system size. For hopping terms whose signs are random (p = 0.5) there 
are such extended domains, since triangles of both phases are distributed all over the 
lattice. In contrast, for small p only small isolated domains are created. 

In this section we address the dependence of the DOS on p. We start by examining 
the dependence of the size of the islands on p. It is natural to expect a percolation 
phase transition, where below a critical probability p c only microscopic islands can be 
found in the lattice. The probability distribution to find an island with area A is then 
exponential with A, where the typical area depends on p and is independent of the 
system size L. Above p c a macroscopic island forms; thus the probability distribution 
is concentrated around the macroscopic value, which of course scales as L 2 . 

This expectation is borne out by our numerical analysis. Figure [8] depicts the 
distribution of the islands area A, normalized to the system's size L 2 , as a function of p 
for an L = 100 lattice. For large p the distribution appears to be concentrated around 
a mean value, which was numerically verified to scale as L 2 . In contrast, for small p the 
distribution is approximately exponential in a manner that was found to be independent 
of L. The transition takes place at p c ~ 0.15. The width of the transition between the 
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microscopic and macroscopic distributions, which takes place here in 0.10 < p < 0.15, 
gets sharper towards p c while increasing L. 

Since flipping a hopping term flips the fluxes of two triangles, islands with even values 
of A are more common than islands with odd values of A. Figure [8] shows that in spite 
of the different distribution of the even and odd values of A, both are exponentially 
distributed below p c . 
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Figure 8. The distribution of the normalized area A/L 2 of the N islands as a 
function of the probability of flipping a hopping term p for a 100 x 100 lattice. (a,b) 
The probability of having an island with area A/L 2 for small p. The probability is 
approximately exponentially distributed, in both even and odd values of A, regardless 
of the system size. This means that for small p in every realization there are many 
microscopic islands, (c) For large p the probability is concentrated around a mean 
macroscopic value, i.e. there is a single macro-island, (d) The probability of having 
an island with area larger than A/L 2 as a function of p. The curve of probability 0.5 
indicates the characteristic largest island. The vertical slope at p c w 0.15 indicates a 
percolation phase transition. The distribution was built upon 100 realizations. 



Figure [8^d) depicts, for every p and normalized area A/L 2 , the probability that there 
are islands larger in area than A/L 2 . The curve for which this probability is 0.5 is 
the typical size of the largest island for a particular disorder realization. We denote 
this curve by A max (p). If we measure the flux in the system with reference to the flux 
of the clean system, then since a triangle with N flux is formed by flipping either one 
or all three hopping terms forming the triangle, the mean flux threading the lattice is 
(0) = 2L 2 [3p(l — p) 2 +p 3 }. Above p c the macroscopic island, if it exists, is supposed to 
capture almost all the N flux. Indeed, A max (p > p c ) is excellently approximated by (0). 
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Thus, above p c almost all the triangles with flipped flux are connected. 

Figure [9] depicts two low-energy disorder-averaged quantities for various lattice sizes 
L. The first is (Ei(p)), the energy of the lowest excitation. The second is the zero- 
energy DOS (-Do)j which we define as the number of states within the energy window 
< E < E sat , divided by D sat , where D sat = (Ei(p — > 0.5)) is the disorder averaged 
energy of the lowest excitation for the p = 0.5 case, extracted from figure |3j In order 
to compare systems with different sizes, we normalize (Ei(p)) by E SS± (L), and (D ) by 
D Bat = (D (p->0.5)). 

For small p we find (Ei) to approach the value of the energy gap, as it should. We 
find (-D ) to be polynomial in p for small p, in accordance with our estimate for the 
probability of finding isolated regions of zero-energy states. Above the percolation phase 
transition the mean minimal energy (Ei) and the zero-energy DOS (D ) are expected 
to saturate to a constant value as a function of p. We do indeed observe a saturation of 
(Ei) and (D ) above p c . 

Surprisingly, figure |9](a) shows that the approach to the saturated value is not 
monotonic. The minimal (E\) appears at p e ~ 0.11 < p c , and we find (Ei(p e )) oc L~ 2 ' 88 , 
as shown in figure |9]^c). Moreover, the maximum (D ) also appear at p e . The divergence 
of the DOS at zero energy seems therefore to be strongest before the percolation 
threshold. 

If the percolation phase transition captures the essence of the low-energy physics, we 
naively expect the low-energy states to be extended above p c , while those below p c , 
including p e , to be localized. For any eigenstate of the Hamiltonian Te = Si Q i7i> 
the localization properties of the state can be characterized by the participation ratio 
(PR), defined by PR = (£\ |aj| 4 ) -1 . The PR of a perfect metal scales linearly with the 
system size L 2 , while the PR of an insulator is independent of L. Figure [9](d) depicts 
the disorder-averaged PR of the lowest excitation as a function of the system size, for 
various values of p. As expected, above p c the PR scales faster than L, which indicates 
an extended state, while well below p c the PR is independent on L. However, at p e the 
PR scales as L. This scaling seems to imply that there is a range of p in which the 
low-energy states are string-like in the sense that they are macroscopic only along one 
dimension. Recalling the observed diverging DOS of the edge states along the narrow 
strips, the divergence of (D ) at p e is understood, even if not the precise exponent. 

5. Summary 

In the previous sections we have seen that the zero-energy peak of the DOS for random 
signs in the nearest neighbor hopping terms may have different sources in two limits of 
the probability p for flipping the sign of individual hopping terms. For small p the peak 
is the consequence of the existence of localized zero-energy states and the tunneling 
between them and other localized states. Therefore it appears on a background of the 
zero DOS of the energy gap, and it is proportional to p 3 . In the opposite limit of 
Pc < P < 0.5, where p c ~ 0.15, macroscopic islands of the time-reversed phase cross the 
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Figure 9. (a) The normalized mean minimal energy (Ei)/E sat as a function of 
the probability p, where E sa ^ — (E±(p — > 0.5)). The percolation phase transition is 
manifested as the saturation to constant values for p > p c , but surprisingly the lowest 
energies appear around p c < p c , which means that they belong to localized states, (b) 
The mean zero-energy DOS (Do) = (D(Q < E < E sat )) as a function of p, normalized 
by its value at p = 0.5, D sat , of the same lattices. The saturation for p > p c is expected, 
but the DOS gets its maximal values at p c < p c , and decays towards zero only in a 
polynomial manner, (c) The mean minimal energy at p c as a function of the system 
size is very well fitted to L" 2 88 . (d) The disorder -averaged participation ratio PR of 
the minimal energy state as a function of the system size L, for various values of p. 
For p > p c the state is extended, since PR scales faster than L (denoted by dashed 
line). For p«p c the states are localized, since PR is independent of L. At p sa p e , 
PR ~ L, which implies a ID string-like state. 



lattice, and the gap is closed due to the low-energy modes along the interfaces of these 
islands with the original phase. Since the dispersion of some of these modes is nonlinear, 
they are expected to create zero-energy peak in the DOS. Our numerical results indicate 
a weak divergence of the DOS at zero energy for p = 0.5, and an even stronger divergence 
somewhat below the percolation transition. In this region the islands seem to be of the 
form of strings, and the mentioned nonlinear dispersion is likely to be the source the 
divergence. 

It was mentioned above that the Hamiltonian of our system is purely imaginary 
and anti-symmetric both to charge conjugation and time-reversal. According to the 
standard classification of disordered Hamiltonians (29] it belongs to class D (30l I3T] . 
This class is known to have three distinct phases: thermal metal, thermal insulator and 
thermal quantum Hall (TQH) insulator [321 133]. In the metal phase, the DOS at zero 
energy diverges logarithmically with the system size [32J. Several numerical studies of 
the properties of this class have been carried out using the Cho-Fisher network model 
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[SUSHI [36], which exhibits all of the phase diagram. It was found that in the TQH phase 
there is a region below the insulator-metal transition where the DOS at zero energy 
diverges, but with a nonuniversal exponent, due to the effects of rare configurations of 
disorder (Griffiths phase) |36j. Recently a model of MFs on a square lattice created 
in a p-wave superconductor by a random potential has been proposed [ID] , and the 
divergence in the metal phase has also been observed numerically there. 

Our model is a new realization to the TQH-thermal metal transition, which may be 
physically realizable. We observe the divergence of the DOS both in the metal regime 
and in the analogue to the Griffiths regime, and the microscopic understanding of our 
system provides insights into the physics of these phenomena. 

Note added: After completing the preparation of this manuscript, we became aware 
of a study of a similar system carried out by C Lohmann, AWW Ludwig and S Trebst 
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Appendix A. Minimal zero mode 

In this appendix we prove that flipping three hopping terms with a common vertex in a 
periodic clean lattice results in two localized states with energy which is either zero or 
exponentially small with the system size. The proof also implies that there is no way 
to create zero modes with flipping less than three signs. 

The hopping signs Sy in the Hamiltonian (jl]) define a matrix 5*, which is skew- 
symmetric S T = —S, and therefore has a well defined Pfaffian. The recursive definition 
of the Pfaffian of an N x N matrix is 

TV 

Pf(A) = £(-l)< ai iPf (Ai), (A.l) 



i=i 




= a Pf (0) = 0, (A.2) 

where is the element i, j of A, and An is the matrix A without the 1 st and i th rows 
and column. A useful property of the Pfaffian [39] is that for any matrix B 

Vi(BAB T ) = det(B) ■ Pf(A). (A.3) 

The gauge transformation 7« — > — 7« flips the signs of the i th row and column of S. This 
transformation can be done by a diagonal matrix B with the elements bjj = (— l) Sj \ 
And since det(B) = —1, the gauge transformation flips the sign of Pf(S'). Moreover, 
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swapping labels of two sites, also flips the sign of Pf(5), since it can be performed by 
B = I(n-2)x(n-2) x o~x, where a x is the Pauli matrix which acts in the space of the two 
swapped sites. 

Systematic gauge transformations reveal equivalence relations between apparently 
different lattices. The simplest example is the gauge transformation that swaps between 
the two sublattices, which is depicted in figure [Alj ^a). In a periodic lattice the number 
of rows is even (which means that the two sublattices are indeed equivalent), an even 
number of sites are gauged, and the Pfaffian is unchanged. Rotations by 60° clockwise 
and counterclockwise can be also produced by gauge transformations, as depicted in 
figures jAjl b) and |Al[ c). In principle, such rotations may multiply the Pfaffian by a 
factor of —1. However, three 60° rotations result in swapped sublattices, which does 
not change the Pfaffian, and hence one rotation must leave the Pfaffian unchanged as 
well. Note that in finite lattices the 60° rotations may have 0(e~ L ) corrections, due to 
change of boundary conditions. 

(C\ o ^ o • ^ • 

AAAA/X 
AAAAAA 
X AAAAA/ 
XAAAA/ 
\AA7y 



/ e 

Figure Al. (a) The gauges transformation that swaps between the two sublattices. 
A red cycle means 7^ — > — 7,-; thus all the arrows connected to the site i are flipped. 
(b,c) Gauge transformation which produce 60° (b) clockwise and (c) anticlockwise 
rotations, (d) An arbitrary site is denoted by 1 and its six neighbors by a-f. 
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All these gauge-dependent properties of the Pfaffian show that it is an unphysical 
quantity. However, since Pf(^4) 2 = det(v4), the Hamiltonian H has zero energy states 
iff Pf(S') = 0. Therefore the way to prove that flipping three neighboring arrows gives 
zero-energy states, would be to prove that the Pfaffian vanishes for such a configuration. 

Let us denote an arbitrary lattice site by 1, and its six neighbors by a,b,..,f, as 



depicted in figure |Al[ d). According to definition (A.l) 

/ 

Pf(S) = ^(-l) i s ll Pf(5 li ). (A.4) 

i=a 

The physical meaning of Pf(Sij) is omitting the sites 1 and i. We look for relations 
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between lattices with vacancies of two neighboring sites. These relations will allow us 



to express five of the terms Pf(Sij) in (A.4) in terms of the sixth, say Pf(5*i a 



Due to the horizontal symmetry of our gauge, it is apparent that |Pf(5i&)| = |Pf(5i a )|, 
but since the labeling of the sites is different, they may differ in signs. If, however, we 
'push' the a th row and column of Su, a — b—1 rows and columns upwards, we get exactly 
Si j0 . This process involves a — b — 1 labels swapping; thus Pf(Sib) = (— l) a_b-1 Pf (Si a ). 
Sic and Su are related by swapping the two sublattices; thus |Pf(jSid)| = |Pf(5i c )|. By 
'pushing' the row and columns, we get again Pi(Sid) = (— l) c_d_1 Pf (Si c ). In a similar 
manner we have Pf(Sif) = (— l) e_ ^Pf (Si e ). Furthermore, Si c (Su) is related to Si a 
by the anticlockwise (clockwise) rotation, hence |Pf(Si e )| ~ |Pf(5i c )| ~ |Pf(£i )|. By 
counting the number of sites been gauges, we get Pf(Si c ) = (— l) a ~ c Pf (Si a ) + 0(e~ L ) 
and Pf(Si e ) = (— l) a-e-1 Pf (Sia) + 0(e~ L ). Substituting all these relations in equation 



flA.4| ) yields 

Pf(S) = (-l) a (s la - si b + s le - s u - s u - Si f )Pf(Si a ) + 0(e~ L ). (A.5) 
In the clean lattice si a = — s lb = s lc = — s ld = — s le = —sif, as can be seen in figure 



Al d). Thus the terms sum-up, and |Pf(5)| = 6|Pf(Si a )|. On the other hand, flipping 
three of the Sji's, gives Pf(S') = 0, i.e. zero-energy modes. Due to the locality of this 
manipulation, the zero modes are localized around site 1. On the other hand, flipping 
only one or two terms will not result in zero modes. 
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